library(glm2)
library(mediation)
library(MASS)
library(dplyr)
library(stargazer)
library(httr)
library(ggplot2)
library(paths)
library(gbm)

set.seed(1234)

setwd(dirname(rstudioapi::getSourceEditorContext()$path))

# Iran 1: Import ----------------------------------------------------------
df_i = read.csv("DataFiles\\Experiment1_Iran.csv")

# Factor certainty and insulting
df_i$certainty = as.factor(df_i$certainty)
df_i$insulting = as.factor(df_i$insulting)

# Iran 1: Mediation Analysis with Dispositional Controls --------------------------

# Subset relevant variables and remove NA values
df_dis = df_i %>% dplyr::select(denial, esca_scaled, MA_scaled, GovTrust, 
                                NewsTrust, IntTrust, NC_scaled, Military.Service, Read.FP,
                                reputation_scaled, certainty, insulting)
df_dis = df_dis[complete.cases(df_dis),]

# Run models with mediators as the DV
m_rep = lm(reputation_scaled ~ denial +  MA_scaled + GovTrust + NewsTrust + IntTrust + NC_scaled + Read.FP, data = df_dis)
m_amb = polr(certainty ~ denial +  MA_scaled + GovTrust + NewsTrust + IntTrust + NC_scaled + Read.FP, data = df_dis, method = "logistic", Hess = TRUE)
m_ins = polr(insulting ~ denial +  MA_scaled + GovTrust + NewsTrust + IntTrust + NC_scaled + Read.FP, data = df_dis, method = "logistic", Hess = TRUE)

# Run models with mediators as the IVs
m2_rep = lm(esca_scaled ~ reputation_scaled + denial +  MA_scaled + GovTrust + NewsTrust + IntTrust + NC_scaled + Read.FP, data = df_dis)
m2_amb = lm(esca_scaled ~ certainty + denial +  MA_scaled + GovTrust + NewsTrust + IntTrust + NC_scaled + Read.FP, data = df_dis)
m2_ins = lm(esca_scaled ~ insulting + denial +  MA_scaled + GovTrust + NewsTrust + IntTrust + NC_scaled + Read.FP, data = df_dis)

# Run mediate function for reputation
med.rep_dis <- mediate(m_rep, m2_rep, treat = "denial", mediator = "reputation_scaled", sims = 500)
summary(med.rep_dis)

# Run mediate function for certainty
med.amb_dis <- mediate(m_amb, m2_amb, treat = "denial", mediator = "certainty", sims = 500, boot = FALSE)
summary(med.amb_dis)

# Run mediate function for insult
med.ins_dis <- mediate(m_ins, m2_ins, treat = "denial", mediator = "insulting", sims = 500, boot = FALSE)
summary(med.ins_dis)

# Iran 1: Mediation Analysis with Demographic Controls ----------------------------

# Subset relevant variables and remove NA values
df_dis = df_i %>% dplyr::select(denial, esca_scaled, certainty, insulting, reputation_scaled, age, male, hhi, white, education, republican, democrat)
df_dis = df_dis[complete.cases(df_dis),]

# Run models with mediators as the DV
m_rep = lm(reputation_scaled ~ denial + age + male + hhi + white + education + republican + democrat, data = df_dis)
m_amb = polr(certainty ~ denial + age + male + hhi + white + education + republican + democrat, data = df_dis, method = "logistic", Hess = TRUE)
m_ins = polr(insulting ~ denial + age + male + hhi + white + education + republican + democrat, data = df_dis, method = "logistic", Hess = TRUE)

# Run models with mediators as the IVs
m2_rep = lm(esca_scaled ~ reputation_scaled + denial + age + male + hhi + white + education + republican + democrat, data = df_dis)
m2_amb = lm(esca_scaled ~ certainty + denial + age + male + hhi + white + education + republican + democrat, data = df_dis)
m2_ins = lm(esca_scaled ~ insulting + denial + age + male + hhi + white + education + republican + democrat, data = df_dis)

# Run mediate function for reputation
med.rep_dem <- mediate(m_rep, m2_rep, treat = "denial", mediator = "reputation_scaled", sims = 500)
summary(med.rep_dem)

# Run mediate function for certainty
med.amb_dem <- mediate(m_amb, m2_amb, treat = "denial", mediator = "certainty", sims = 500, boot = FALSE)
summary(med.amb_dem)

# Run mediate function for insult
med.ins_dem <- mediate(m_ins, m2_ins, treat = "denial", mediator = "insulting", sims = 500, boot = FALSE)
summary(med.ins_dem)

# Iran 1: Mediation Analysis Results Table (Table 1 and Table A6) ----------------------------------------

# Create HTML table of full mediation results
med_table = paste('<!DOCTYPE html>
  <html>
  <style>
  table
  th{text-align: left;}
</style>
  <body>
  
  <h2>Iran 1: Mediation Effects</h2>
  
  <table style="width:35%">
  <tr><th colspan="5" style = "border-bottom: 1px solid black"></th></tr>
  <tr align = "left">
    <th style="width:15%"></th>
    <th style="width:15%">Mediator</th>
    <th style="width:15%">ACME</th>
    <th style="width:15%">Direct Effect</th>
    <th style="width:15%">Total Effect</th>
  </tr>
  <tr><th colspan="5" style = "border-bottom: 1px solid black"></th></tr>
  <tr>
    <th rowspan ="3">Dispositional Controls</th> 
    <td >Reputation</td>
    <td>', round(med.rep_dis$d0,2), '(',round(med.rep_dis$d0.ci[1],2),',',round(med.rep_dis$d0.ci[2],2),')','</td>
    <td>', round(med.rep_dis$z0,2),'(',round(med.rep_dis$z0.ci[1],2),',',round(med.rep_dis$z0.ci[2],2),')','</td>
    <td>', round(med.rep_dis$tau.coef,2),'(',round(med.rep_dis$tau.ci[1],2),',',round(med.rep_dis$tau.ci[2],2),')','</td>
  </tr>
  <tr style="background-color: #D6EEEE">
    <td>Certainty</td>
    <td>', round(med.amb_dis$d0,2), '(',round(med.amb_dis$d0.ci[1],2),',',round(med.amb_dis$d0.ci[2],2),')','</td>
    <td>', round(med.amb_dis$z0,2),'(',round(med.amb_dis$z0.ci[1],2),',',round(med.amb_dis$z0.ci[2],2),')','</td>
    <td>', round(med.amb_dis$tau.coef,2),'(',round(med.amb_dis$tau.ci[1],2),',',round(med.amb_dis$tau.ci[2],2),')','</td>
  </tr>
  <tr>
    <td>Insulting</td>
   <td>', round(med.ins_dis$d0,2),'(',round(med.ins_dis$d0.ci[1],2),',',round(med.ins_dis$d0.ci[2],2),')','</td>
    <td>', round(med.ins_dis$z0,2),'(',round(med.ins_dis$z0.ci[1],2),',',round(med.ins_dis$z0.ci[2],2),')','</td>
    <td>', round(med.ins_dis$tau.coef,2),'(',round(med.ins_dis$tau.ci[1],2),',',round(med.ins_dis$tau.ci[2],2),')','</td>
  </tr>
  <tr><th colspan="6" style = "border-bottom: 1px solid black"></th></tr>
  <tr>
    <th rowspan ="3">Demographic Controls</th> 
    <td>Reputation</td>
    <td>', round(med.rep_dem$d0,2), '(',round(med.rep_dem$d0.ci[1],2),',',round(med.rep_dem$d0.ci[2],2),')','</td>
    <td>', round(med.rep_dem$z0,2),'(',round(med.rep_dem$z0.ci[1],2),',',round(med.rep_dem$z0.ci[2],2),')','</td>
    <td>', round(med.rep_dem$tau.coef,2),'(',round(med.rep_dem$tau.ci[1],2),',',round(med.rep_dem$tau.ci[2],2),')','</td>
  </tr>
  <tr style="background-color: #D6EEEE">
    <td>Certainty</td>
    <td>', round(med.amb_dem$d0,2), '(',round(med.amb_dem$d0.ci[1],2),',',round(med.amb_dem$d0.ci[2],2),')','</td>
    <td>', round(med.amb_dem$z0,2),'(',round(med.amb_dem$z0.ci[1],2),',',round(med.amb_dem$z0.ci[2],2),')','</td>
    <td>', round(med.amb_dem$tau.coef,2),'(',round(med.amb_dem$tau.ci[1],2),',',round(med.amb_dem$tau.ci[2],2),')','</td>
  </tr>
  <tr>
    <td>Insulting</td>
    <td>', round(med.ins_dem$d0,2), '(',round(med.ins_dem$d0.ci[1],2),',',round(med.ins_dem$d0.ci[2],2),')','</td>
    <td>', round(med.ins_dem$z0,2),'(',round(med.ins_dem$z0.ci[1],2),',',round(med.ins_dem$z0.ci[2],2),')','</td>
    <td>', round(med.ins_dem$tau.coef,2),'(',round(med.ins_dem$tau.ci[1],2),',',round(med.ins_dem$tau.ci[2],2),')','</td>
  </tr>
  <tr><th colspan="5" style = "border-bottom: 1px solid black"></th></tr>
  </table>
  
  <p></p>
  
  </body>
  </html>', sep =" ")
write(med_table, "Tables&Figures/Iran 1 Mediation Results (Tables 1 and A6).html")
BROWSE("Tables&Figures/Iran 1 Mediation Results (Tables 1 and A6).html")

# We calculate the proportion mediated manually by dividing the ACME by the Total Effect.

# Iran 1: Mediation Analysis Controlling for Causal Dependency (Table A7) --------------------

# Multi-mediation Imai and Yamamoto
# certainty
df_dis = df_i %>% dplyr::select(denial, esca_scaled, MA_scaled, GovTrust, 
                                NewsTrust, IntTrust, NC_scaled, Read.FP,
                                reputation_scaled, certainty, insulting, age, male, 
                                hhi, white, education, republican, democrat)
df_dis = df_dis[complete.cases(df_dis),]


df_dis$certainty = as.numeric(df_dis$certainty)
df_dis$insulting = as.numeric(df_dis$insulting)

covariates = c('MA_scaled', 'GovTrust', 'NewsTrust', 'IntTrust', 'NC_scaled', 'Read.FP')
meds = c('reputation_scaled','insulting')
mm_amb <- multimed("esca_scaled", "certainty", meds, "denial", covariates, 
                   data=df_dis, sims=1000)
summary(mm_amb)

# Reputation
covariates = c('MA_scaled', 'GovTrust', 'NewsTrust', 'IntTrust', 'NC_scaled', 'Read.FP')
meds = c('certainty','insulting')
mm_rep <- multimed("esca_scaled", "reputation_scaled", meds, "denial", covariates, 
                   data=df_dis, sims=1000)
summary(mm_rep)

# Insulting
covariates = c('MA_scaled', 'GovTrust', 'NewsTrust', 'IntTrust', 'NC_scaled', 'Read.FP')
meds = c('certainty', 'reputation_scaled')
mm_ins <- multimed("esca_scaled", "insulting", meds, "denial", covariates, 
                   data=df_dis, sims=1000)
summary(mm_ins)


# Qatar 2: Import ---------------------------------------------------------
df_i = read.csv("DataFiles\\Experiment2_Qatar.csv")

# Factor certainty and insulting
df_i$certainty = as.factor(df_i$certainty)
df_i$insulting = as.factor(df_i$insulting)

# Qatar 2: Mediation Analysis with Dispositional Controls --------------------------

# Subset relevant variables and remove NA values
df_dis = df_i %>% dplyr::select(denial, esca_scaled, MA_scaled, GovTrust, 
                                NewsTrust, IntTrust, NC_scaled, Military.Service, Read.FP,
                                reputation_scaled, certainty, insulting)
df_dis = df_dis[complete.cases(df_dis),]

# Run models with mediators as the DV
m_rep = lm(reputation_scaled ~ denial +  MA_scaled + GovTrust + NewsTrust + IntTrust + NC_scaled + Read.FP, data = df_dis)
m_amb = polr(certainty ~ denial +  MA_scaled + GovTrust + NewsTrust + IntTrust + NC_scaled + Read.FP, data = df_dis, method = "logistic", Hess = TRUE)
m_ins = polr(insulting ~ denial +  MA_scaled + GovTrust + NewsTrust + IntTrust + NC_scaled + Read.FP, data = df_dis, method = "logistic", Hess = TRUE)

# Run models with mediators as the IVs
m2_rep = lm(esca_scaled ~ reputation_scaled + denial +  MA_scaled + GovTrust + NewsTrust + IntTrust + NC_scaled + Read.FP, data = df_dis)
m2_amb = lm(esca_scaled ~ certainty + denial +  MA_scaled + GovTrust + NewsTrust + IntTrust + NC_scaled + Read.FP, data = df_dis)
m2_ins = lm(esca_scaled ~ insulting + denial +  MA_scaled + GovTrust + NewsTrust + IntTrust + NC_scaled + Read.FP, data = df_dis)

# Run mediate function for reputation
med.rep_dis <- mediate(m_rep, m2_rep, treat = "denial", mediator = "reputation_scaled", robustSE = TRUE, sims = 500)
summary(med.rep_dis)

# Run mediate function for certainty
med.amb_dis <- mediate(m_amb, m2_amb, treat = "denial", mediator = "certainty", sims = 500, boot = FALSE)
summary(med.amb_dis)

# Run mediate function for insult
med.ins_dis <- mediate(m_ins, m2_ins, treat = "denial", mediator = "insulting", sims = 500, boot = FALSE)
summary(med.ins_dis)

# Qatar 2: Mediation Analysis with Demographic Controls ----------------------------

# Subset relevant variables and remove NA values
df_dis = df_i %>% dplyr::select(denial, esca_scaled, certainty, insulting, reputation_scaled, age, male, hhi, white, education, republican, democrat)
df_dis = df_dis[complete.cases(df_dis),]

# Run models with mediators as the DV
m_rep = lm(reputation_scaled ~ denial + age + male + hhi + white + education + republican + democrat, data = df_dis)
m_amb = polr(certainty ~ denial + age + male + hhi + white + education + republican + democrat, data = df_dis, method = "logistic", Hess = TRUE)
m_ins = polr(insulting ~ denial + age + male + hhi + white + education + republican + democrat, data = df_dis, method = "logistic", Hess = TRUE)

# Run models with mediators as the IVs
m2_rep = lm(esca_scaled ~ reputation_scaled + denial + age + male + hhi + white + education + republican + democrat, data = df_dis)
m2_amb = lm(esca_scaled ~ certainty + denial + age + male + hhi + white + education + republican + democrat, data = df_dis)
m2_ins = lm(esca_scaled ~ insulting + denial + age + male + hhi + white + education + republican + democrat, data = df_dis)

# Run mediate function for reputation
med.rep_dem <- mediate(m_rep, m2_rep, treat = "denial", mediator = "reputation_scaled", robustSE = TRUE, sims = 500)
summary(med.rep_dem)

# Run mediate function for certainty
med.amb_dem <- mediate(m_amb, m2_amb, treat = "denial", mediator = "certainty", sims = 500, boot = FALSE)
summary(med.amb_dem)

# Run mediate function for insult
med.ins_dem <- mediate(m_ins, m2_ins, treat = "denial", mediator = "insulting", sims = 500, boot = FALSE)
summary(med.ins_dem)

# Qatar 2: Mediation Analysis Results Table (Table 2 and Table A13) ----------------------------------------

# Create HTML table of full mediation results
med_table = paste('<!DOCTYPE html>
  <html>
  <style>
  table
  th{text-align: left;}
</style>
  <body>
  
  <h2>Qatar 2: Mediation Effects</h2>
  
  <table style="width:35%">
  <tr><th colspan="5" style = "border-bottom: 1px solid black"></th></tr>
  <tr align = "left">
    <th style="width:15%"></th>
    <th style="width:15%">Mediator</th>
    <th style="width:15%">ACME</th>
    <th style="width:15%">Direct Effect</th>
    <th style="width:15%">Total Effect</th>
  </tr>
  <tr><th colspan="5" style = "border-bottom: 1px solid black"></th></tr>
  <tr>
    <th rowspan ="3">Dispositional Controls</th> 
    <td >Reputation</td>
    <td>', round(med.rep_dis$d0,2), '(',round(med.rep_dis$d0.ci[1],2),',',round(med.rep_dis$d0.ci[2],2),')','</td>
    <td>', round(med.rep_dis$z0,2),'(',round(med.rep_dis$z0.ci[1],2),',',round(med.rep_dis$z0.ci[2],2),')','</td>
    <td>', round(med.rep_dis$tau.coef,2),'(',round(med.rep_dis$tau.ci[1],2),',',round(med.rep_dis$tau.ci[2],2),')','</td>
  </tr>
  <tr style="background-color: #D6EEEE">
    <td>Certainty</td>
    <td>', round(med.amb_dis$d0,2), '(',round(med.amb_dis$d0.ci[1],2),',',round(med.amb_dis$d0.ci[2],2),')','</td>
    <td>', round(med.amb_dis$z0,2),'(',round(med.amb_dis$z0.ci[1],2),',',round(med.amb_dis$z0.ci[2],2),')','</td>
    <td>', round(med.amb_dis$tau.coef,2),'(',round(med.amb_dis$tau.ci[1],2),',',round(med.amb_dis$tau.ci[2],2),')','</td>
  </tr>
  <tr>
    <td>Insulting</td>
   <td>', round(med.ins_dis$d0,2),'(',round(med.ins_dis$d0.ci[1],2),',',round(med.ins_dis$d0.ci[2],2),')','</td>
    <td>', round(med.ins_dis$z0,2),'(',round(med.ins_dis$z0.ci[1],2),',',round(med.ins_dis$z0.ci[2],2),')','</td>
    <td>', round(med.ins_dis$tau.coef,2),'(',round(med.ins_dis$tau.ci[1],2),',',round(med.ins_dis$tau.ci[2],2),')','</td>
  </tr>
  <tr><th colspan="6" style = "border-bottom: 1px solid black"></th></tr>
  <tr>
    <th rowspan ="3">Demographic Controls</th> 
    <td>Reputation</td>
    <td>', round(med.rep_dem$d0,2), '(',round(med.rep_dem$d0.ci[1],2),',',round(med.rep_dem$d0.ci[2],2),')','</td>
    <td>', round(med.rep_dem$z0,2),'(',round(med.rep_dem$z0.ci[1],2),',',round(med.rep_dem$z0.ci[2],2),')','</td>
    <td>', round(med.rep_dem$tau.coef,2),'(',round(med.rep_dem$tau.ci[1],2),',',round(med.rep_dem$tau.ci[2],2),')','</td>
  </tr>
  <tr style="background-color: #D6EEEE">
    <td>Certainty</td>
    <td>', round(med.amb_dem$d0,2), '(',round(med.amb_dem$d0.ci[1],2),',',round(med.amb_dem$d0.ci[2],2),')','</td>
    <td>', round(med.amb_dem$z0,2),'(',round(med.amb_dem$z0.ci[1],2),',',round(med.amb_dem$z0.ci[2],2),')','</td>
    <td>', round(med.amb_dem$tau.coef,2),'(',round(med.amb_dem$tau.ci[1],2),',',round(med.amb_dem$tau.ci[2],2),')','</td>
  </tr>
  <tr>
    <td>Insulting</td>
    <td>', round(med.ins_dem$d0,2), '(',round(med.ins_dem$d0.ci[1],2),',',round(med.ins_dem$d0.ci[2],2),')','</td>
    <td>', round(med.ins_dem$z0,2),'(',round(med.ins_dem$z0.ci[1],2),',',round(med.ins_dem$z0.ci[2],2),')','</td>
    <td>', round(med.ins_dem$tau.coef,2),'(',round(med.ins_dem$tau.ci[1],2),',',round(med.ins_dem$tau.ci[2],2),')','</td>
  </tr>
  <tr><th colspan="5" style = "border-bottom: 1px solid black"></th></tr>
  </table>
  
  <p></p>
  
  </body>
  </html>', sep =" ")
write(med_table, "Tables&Figures/Qatar 2 Mediation Results (Tables 2 and A13).html")
BROWSE("Tables&Figures/Qatar 2 Mediation Results (Tables 2 and A13).html")

# We calculate the proportion mediated manually by dividing the ACME by the Total Effect.

# Qatar 2: Mediation Analysis Controlling for Causal Dependency (Table A14) --------------------

# Multi-mediation Imai and Yamamoto
# certainty
df_dis = df_i %>% dplyr::select(denial, esca_scaled, MA_scaled, GovTrust, 
                                NewsTrust, IntTrust, NC_scaled, Read.FP,
                                reputation_scaled, certainty, insulting, age, male, 
                                hhi, white, education, republican, democrat)
df_dis = df_dis[complete.cases(df_dis),]


df_dis$certainty = as.numeric(df_dis$certainty)
df_dis$insulting = as.numeric(df_dis$insulting)

covariates = c('MA_scaled', 'GovTrust', 'NewsTrust', 'IntTrust', 'NC_scaled', 'Read.FP')
meds = c('reputation_scaled','insulting')
mm_amb <- multimed("esca_scaled", "certainty", meds, "denial", covariates, 
                   data=df_dis, sims=1000)
summary(mm_amb)

# Reputation
covariates = c('MA_scaled', 'GovTrust', 'NewsTrust', 'IntTrust', 'NC_scaled', 'Read.FP')
meds = c('certainty','insulting')
mm_rep <- multimed("esca_scaled", "reputation_scaled", meds, "denial", covariates, 
                   data=df_dis, sims=1000)
summary(mm_rep)

# Insulting
covariates = c('MA_scaled', 'GovTrust', 'NewsTrust', 'IntTrust', 'NC_scaled', 'Read.FP')
meds = c('certainty', 'reputation_scaled')
mm_ins <- multimed("esca_scaled", "insulting", meds, "denial", covariates, 
                   data=df_dis, sims=1000)
summary(mm_ins)




# Iran 3: Import ----------------------------------------------------------
df_i = read.csv("DataFiles\\Experiment3_Iran.csv")

# Factor certainty and insulting
df_i$certainty = as.factor(df_i$certainty)
df_i$insulting = as.factor(df_i$insulting)

# Iran 3: Mediation Analysis with Dispositional Controls --------------------------

# Subset relevant variables and remove NA values
df_dis = df_i %>% dplyr::select(denial, esca_scaled, MA_scaled, GovTrust, 
                                NewsTrust, IntTrust, NC_scaled, Read.FP,
                                reputation_scaled, certainty, insulting, age, male, 
                                hhi, white, education, republican, democrat)
df_dis = df_dis[complete.cases(df_dis),]

# Run models with mediators as the DV
m_rep = lm(reputation_scaled ~ denial +  MA_scaled + GovTrust + NewsTrust + IntTrust + NC_scaled + Read.FP, data = df_dis)
m_amb = polr(certainty ~ denial +  MA_scaled + GovTrust + NewsTrust + IntTrust + NC_scaled + Read.FP, data = df_dis, method = "logistic", Hess = TRUE)
m_ins = polr(insulting ~ denial +  MA_scaled + GovTrust + NewsTrust + IntTrust + NC_scaled + Read.FP, data = df_dis, method = "logistic", Hess = TRUE)

# Run models with mediators as the IVs
m2_rep = lm(esca_scaled ~ reputation_scaled + denial +  MA_scaled + GovTrust + NewsTrust + IntTrust + NC_scaled + Read.FP, data = df_dis)
m2_amb = lm(esca_scaled ~ certainty + denial +  MA_scaled + GovTrust + NewsTrust + IntTrust + NC_scaled + Read.FP, data = df_dis)
m2_ins = lm(esca_scaled ~ insulting + denial +  MA_scaled + GovTrust + NewsTrust + IntTrust + NC_scaled + Read.FP, data = df_dis)

# Run mediate function for reputation
med.rep_dis <- mediate(m_rep, m2_rep, treat = "denial", mediator = "reputation_scaled", sims = 500, boot = FALSE)
summary(med.rep_dis)

# Run mediate function for certainty
med.amb_dis <- mediate(m_amb, m2_amb, treat = "denial", mediator = "certainty", sims = 500, boot = FALSE)
summary(med.amb_dis)

# Run mediate function for insult
med.ins_dis <- mediate(m_ins, m2_ins, treat = "denial", mediator = "insulting", sims = 500, boot = FALSE)
summary(med.ins_dis)

# Iran 3: Mediation Analysis with Demographic Controls ----------------------------

# Subset relevant variables and remove NA values
df_dis = df_i %>% dplyr::select(denial, esca_scaled, certainty, insulting, reputation_scaled, age, male, hhi, white, education, republican, democrat)
df_dis = df_dis[complete.cases(df_dis),]

# Run models with mediators as the DV
m_rep = lm(reputation_scaled ~ denial + age + male + hhi + white + education + republican + democrat, data = df_dis)
m_amb = polr(certainty ~ denial + age + male + hhi + white + education + republican + democrat, data = df_dis, method = "logistic", Hess = TRUE)
m_ins = polr(insulting ~ denial + age + male + hhi + white + education + republican + democrat, data = df_dis, method = "logistic", Hess = TRUE)

# Run models with mediators as the IVs
m2_rep = lm(esca_scaled ~ reputation_scaled + denial + age + male + hhi + white + education + republican + democrat, data = df_dis)
m2_amb = lm(esca_scaled ~ certainty + denial + age + male + hhi + white + education + republican + democrat, data = df_dis)
m2_ins = lm(esca_scaled ~ insulting + denial + age + male + hhi + white + education + republican + democrat, data = df_dis)

# Run mediate function for reputation
med.rep_dem <- mediate(m_rep, m2_rep, treat = "denial", mediator = "reputation_scaled", sims = 500)
summary(med.rep_dem)

# Run mediate function for certainty
med.amb_dem <- mediate(m_amb, m2_amb, treat = "denial", mediator = "certainty", sims = 500, boot = FALSE)
summary(med.amb_dem)

# Run mediate function for insult
med.ins_dem <- mediate(m_ins, m2_ins, treat = "denial", mediator = "insulting", sims = 500, boot = FALSE)
summary(med.ins_dem)

# Iran 3: Mediation Analysis Results Table (Table 3 and Table A20) ----------------------------------------

# Create HTML table of full mediation results
med_table = paste('<!DOCTYPE html>
  <html>
  <style>
  table
  th{text-align: left;}
</style>
  <body>
  
  <h2>Iran 3: Mediation Effects</h2>
  
  <table style="width:35%">
  <tr><th colspan="5" style = "border-bottom: 1px solid black"></th></tr>
  <tr align = "left">
    <th style="width:15%"></th>
    <th style="width:15%">Mediator</th>
    <th style="width:15%">ACME</th>
    <th style="width:15%">Direct Effect</th>
    <th style="width:15%">Total Effect</th>
  </tr>
  <tr><th colspan="5" style = "border-bottom: 1px solid black"></th></tr>
  <tr>
    <th rowspan ="3">Dispositional Controls</th> 
    <td >Reputation</td>
    <td>', round(med.rep_dis$d0,2), '(',round(med.rep_dis$d0.ci[1],2),',',round(med.rep_dis$d0.ci[2],2),')','</td>
    <td>', round(med.rep_dis$z0,2),'(',round(med.rep_dis$z0.ci[1],2),',',round(med.rep_dis$z0.ci[2],2),')','</td>
    <td>', round(med.rep_dis$tau.coef,2),'(',round(med.rep_dis$tau.ci[1],2),',',round(med.rep_dis$tau.ci[2],2),')','</td>
  </tr>
  <tr style="background-color: #D6EEEE">
    <td>Certainty</td>
    <td>', round(med.amb_dis$d0,2), '(',round(med.amb_dis$d0.ci[1],2),',',round(med.amb_dis$d0.ci[2],2),')','</td>
    <td>', round(med.amb_dis$z0,2),'(',round(med.amb_dis$z0.ci[1],2),',',round(med.amb_dis$z0.ci[2],2),')','</td>
    <td>', round(med.amb_dis$tau.coef,2),'(',round(med.amb_dis$tau.ci[1],2),',',round(med.amb_dis$tau.ci[2],2),')','</td>
  </tr>
  <tr>
    <td>Insulting</td>
   <td>', round(med.ins_dis$d0,2),'(',round(med.ins_dis$d0.ci[1],2),',',round(med.ins_dis$d0.ci[2],2),')','</td>
    <td>', round(med.ins_dis$z0,2),'(',round(med.ins_dis$z0.ci[1],2),',',round(med.ins_dis$z0.ci[2],2),')','</td>
    <td>', round(med.ins_dis$tau.coef,2),'(',round(med.ins_dis$tau.ci[1],2),',',round(med.ins_dis$tau.ci[2],2),')','</td>
  </tr>
  <tr><th colspan="6" style = "border-bottom: 1px solid black"></th></tr>
  <tr>
    <th rowspan ="3">Demographic Controls</th> 
    <td>Reputation</td>
    <td>', round(med.rep_dem$d0,2), '(',round(med.rep_dem$d0.ci[1],2),',',round(med.rep_dem$d0.ci[2],2),')','</td>
    <td>', round(med.rep_dem$z0,2),'(',round(med.rep_dem$z0.ci[1],2),',',round(med.rep_dem$z0.ci[2],2),')','</td>
    <td>', round(med.rep_dem$tau.coef,2),'(',round(med.rep_dem$tau.ci[1],2),',',round(med.rep_dem$tau.ci[2],2),')','</td>
  </tr>
  <tr style="background-color: #D6EEEE">
    <td>Certainty</td>
    <td>', round(med.amb_dem$d0,2), '(',round(med.amb_dem$d0.ci[1],2),',',round(med.amb_dem$d0.ci[2],2),')','</td>
    <td>', round(med.amb_dem$z0,2),'(',round(med.amb_dem$z0.ci[1],2),',',round(med.amb_dem$z0.ci[2],2),')','</td>
    <td>', round(med.amb_dem$tau.coef,2),'(',round(med.amb_dem$tau.ci[1],2),',',round(med.amb_dem$tau.ci[2],2),')','</td>
  </tr>
  <tr>
    <td>Insulting</td>
    <td>', round(med.ins_dem$d0,2), '(',round(med.ins_dem$d0.ci[1],2),',',round(med.ins_dem$d0.ci[2],2),')','</td>
    <td>', round(med.ins_dem$z0,2),'(',round(med.ins_dem$z0.ci[1],2),',',round(med.ins_dem$z0.ci[2],2),')','</td>
    <td>', round(med.ins_dem$tau.coef,2),'(',round(med.ins_dem$tau.ci[1],2),',',round(med.ins_dem$tau.ci[2],2),')','</td>
  </tr>
  <tr><th colspan="5" style = "border-bottom: 1px solid black"></th></tr>
  </table>
  
  <p></p>
  
  </body>
  </html>', sep =" ")
write(med_table, "Tables&Figures/Iran 3 Mediation Results (Tables 3 and A20).html")
BROWSE("Tables&Figures/Iran 3 Mediation Results (Tables 3 and A20).html")

# We calculate the proportion mediated manually by dividing the ACME by the Total Effect, except for
# certainty. Because the ACME for certainty is larger than the total effect, we calculate the proportion
# mediated by dividing the total effect by the ACME. See the footnote in the manuscript.

# Iran 3: Mediation Analysis Controlling for Causal Dependency (Table A21) --------------------

# Multi-mediation Imai and Yamamoto
# certainty
df_dis = df_i %>% dplyr::select(denial, esca_scaled, MA_scaled, GovTrust, 
                                NewsTrust, IntTrust, NC_scaled, Read.FP,
                                reputation_scaled, certainty, insulting, age, male, 
                                hhi, white, education, republican, democrat)
df_dis = df_dis[complete.cases(df_dis),]

df_dis$certainty = as.numeric(df_dis$certainty)
df_dis$insulting = as.numeric(df_dis$insulting)

covariates = c('MA_scaled', 'GovTrust', 'NewsTrust', 'IntTrust', 'NC_scaled', 'Read.FP')
meds = c('reputation_scaled','insulting')
mm_amb <- multimed("esca_scaled", "certainty", meds, "denial", covariates, 
                   data=df_dis, sims=1000)
summary(mm_amb)

# Reputation
covariates = c('MA_scaled', 'GovTrust', 'NewsTrust', 'IntTrust', 'NC_scaled', 'Read.FP')
meds = c('certainty','insulting')
mm_rep <- multimed("esca_scaled", "reputation_scaled", meds, "denial", covariates, 
                   data=df_dis, sims=1000)
summary(mm_rep)

# Insulting
covariates = c('MA_scaled', 'GovTrust', 'NewsTrust', 'IntTrust', 'NC_scaled', 'Read.FP')
meds = c('certainty', 'reputation_scaled')
mm_ins <- multimed("esca_scaled", "insulting", meds, "denial", covariates, 
                   data=df_dis, sims=1000)
summary(mm_ins)


